Setup

library(sp)
library(raster)
library(R.matlab)
library(plyr)
library(data.table)
library(maptools)
library(rgdal)
library(spatstat)

Data Proccessing

setwd("~/Desktop/Fall_2018_Classes/GIS_713/Final_Project")
library(plyr)
library(data.table)
CHIPS_df = read.table("chips.csv", header = TRUE, row.names=NULL, sep=",")
CHIPS_df <- CHIPS_df[!CHIPS_df$lat < 500000.00, ]
CHIPS_df <- CHIPS_df[!CHIPS_df$long < 50000.00, ]
# correcting odd lat coords
CHIPS_df$lat <- CHIPS_df$lat / 10000
CHIPS_df$long <- CHIPS_df$long / 10000
# And a function that shifts vectors conviniently (we gonna need this later):
shift.vec <- function (vec, shift) {
  if(length(vec) <= abs(shift)) {
    rep(NA ,length(vec))
  }else{
    if (shift >= 0) {
      c(rep(NA, shift), vec[1:(length(vec)-shift)]) }
    else {
      c(vec[(abs(shift)+1):length(vec)], rep(NA, abs(shift))) } }
  }
# Now, let’s calculate the distances between successive positions and the respective speed in this segment.
# Shift vectors for lat and lon so that each row also contains the next position.
CHIPS_df$lat.p1 <- shift.vec(CHIPS_df$lat, -1)
CHIPS_df$lon.p1 <- shift.vec(CHIPS_df$long, -1)
# Calculate distances (in metres) using the function pointDistance from the ‘raster’ package.
# Parameter ‘lonlat’ has to be TRUE!
CHIPS_df$dist.to.prev <- apply(CHIPS_df, 1, FUN = function (row) {
  pointDistance(c(as.numeric(row["lat.p1"]), as.numeric(row["long.p1"])), 
                c(as.numeric(row["lat"]), as.numeric(row["long"])), lonlat = T)
})
# Transform the column ‘time’ so that R knows how to interpret it.
CHIPS_df$time_new <- strptime(CHIPS_df$initial_time_stamp_mat, format="%m/%d/%Y %H:%M")
# Shift the time vector, too.
CHIPS_df$time.p1 <- shift.vec(CHIPS_df$time_new, -1)
# Calculate the number of seconds between two positions.
CHIPS_df$time.diff.to.prev <- as.numeric(difftime(CHIPS_df$time.p1, CHIPS_df$time_new))
# Calculate metres per seconds, kilometres per hour and two LOWESS smoothers to get rid of some noise.
CHIPS_df$speed.m.per.sec <- CHIPS_df$dist.to.prev / CHIPS_df$time.diff.to.prev
CHIPS_df$speed.km.per.h <- CHIPS_df$speed.m.per.sec * 3.6
CHIPS_df$speed.km.per.h <- ifelse(is.na(CHIPS_df$speed.km.per.h), 0, CHIPS_df$speed.km.per.h)
CHIPS_df$lowess.speed <- lowess(CHIPS_df$speed2, f = 0.2)$y
CHIPS_df$lowess.alt <- lowess(CHIPS_df$altitude, f = 0.2)$y
CHIPS_df$lowess.conduct <- lowess(CHIPS_df$conductance_z, f = 0.2)$y

setnames(CHIPS_df, "long", "lon")
pt1 <- CHIPS_df[CHIPS_df$participant == 1, ]
pt2 <- CHIPS_df[CHIPS_df$participant == 2, ]
pt3 <- CHIPS_df[CHIPS_df$participant == 3, ]
pt4 <- CHIPS_df[CHIPS_df$participant == 4, ]
pt5 <- CHIPS_df[CHIPS_df$participant == 5, ]
pt6 <- CHIPS_df[CHIPS_df$participant == 6, ]
pt7 <- CHIPS_df[CHIPS_df$participant == 7, ]
pt8 <- CHIPS_df[CHIPS_df$participant == 8, ]
pt9 <- CHIPS_df[CHIPS_df$participant == 9, ]
pt10 <- CHIPS_df[CHIPS_df$participant == 10, ]
pt11 <- CHIPS_df[CHIPS_df$participant == 11, ]

Exploratory Data Visualizations

GPS

# Now, let’s plot all the stuff!
# Plot elevations and smoother
layout(matrix(1:3, nrow=3))
plot(CHIPS_df$altitude, type = "l", bty = "n", xaxt = "n", lwd= 3, 
     ylab = "Elevation", 
     xlab = "", col = "grey60")
lines(CHIPS_df$lowess.alt, col = "green", lwd = 3)
abline(h = mean(CHIPS_df$altitude), lty = 2, lwd = 3, col = "green")
legend(x="bottomright", legend = c("GPS elevation", "LOWESS elevation", 
                                "Mean elevation"),
       col = c("grey60", "green", "green"), lwd = c(2,4,2), lty = c(1,2,2),
       bty = "n")
# Plot speeds and smoother
plot(CHIPS_df$speed2, type = "l", bty = "n", lwd= 3, xaxt = "n", 
     ylab = "Speed (km/h)", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.speed, col = "red", lwd = 3)
abline(h = mean(CHIPS_df$speed2), lty = 2, lwd = 3, col = "red")
legend(x="topright", legend = c("GPS speed", "LOWESS speed", 
                                   "Mean speed"),
       col = c("grey60", "red", "red"), lwd = c(2,4,2), lty = c(1,2,2), 
       bty = "n")
# Plot conductnace and smoother
plot(CHIPS_df$conductance_z, type = "l", bty = "n", lwd= 3, xaxt = "n", 
     ylab = "Skin Conductance", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.conduct, col = "blue", lwd = 3)
abline(h = mean(CHIPS_df$conductance_z), lty = 2, lwd = 3,  col = "blue")
legend(x="topright",
       legend = c("Conductance", "LOWESS conductance", "Mean conductance"),
       col = c("grey60", "blue", "blue"), lwd = c(2,4,2), lty = c(1,2,2),
       bty = "n")

par(mar=c(5, 4, 4, 2) + 0.1)

Skin Conductance

# Plotting the elevation and timestamp of each waypoint using ggplot, allows us to visualise the hike towards the valley of Grauson. As a preparatory step, we use the ymd_hms() function from the lubridate library to convert the string representating the timestamp to a proper R time-object. As to not confuse ggplot, we also do not pass the SpatialPointsDataFrame-object directly, but convert it to a regular dataframe with as.data.frame():
if(!require(lubridate)) install_github("rstudio/lubridate")
if(!require(ggplot2)) install_github("rstudio/ggplot2")
if(!require(gridExtra)) install_github("rstudio/gridExtra")
# plot of time and elevation, colored by skin conductance
time_ele_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
            aes(x=time, y=altitude, color = conductance_z)) +
            scale_color_gradient(low="blue", high="red") +
            geom_point(alpha = 0.01, size = 4) + 
            labs(x='Cycling time', y='Elevations (meters)')
# plot of time and speed, colored by skin conductance
time_speed_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
            aes(x=time, y=speed2, color = conductance_z)) +
            scale_color_gradient(low="blue", high="red") +
            geom_point(alpha = 0.08, size = 4) + 
  labs(x='Cycling time', y='Speed (km/h)')
grid.arrange(time_ele_conduct_plot, time_speed_conduct_plot, nrow=2)

Spatial Data Processing

shp_dsn <- "~/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles"
landcover <- readOGR(path.expand(shp_dsn), 'NL006L3_TILBURG_UA2012')
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/garrettmillar/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles", layer: "NL006L3_TILBURG_UA2012"
## with 12086 features
## It has 11 fields
# Projection
landcover <- spTransform(landcover, CRS('+proj=longlat'))
# Conversion into SpatialPoints
coordinates(CHIPS_df) <- ~ lon + lat
coordinates(pt1) <- ~ lon + lat
coordinates(pt2) <- ~ lon + lat
coordinates(pt3) <- ~ lon + lat
coordinates(pt4) <- ~ lon + lat
coordinates(pt5) <- ~ lon + lat
coordinates(pt6) <- ~ lon + lat
coordinates(pt7) <- ~ lon + lat
coordinates(pt8) <- ~ lon + lat
coordinates(pt9) <- ~ lon + lat
coordinates(pt10) <- ~ lon + lat
coordinates(pt11) <- ~ lon + lat
# Setting default projection
proj4string(CHIPS_df) <- CRS('+proj=longlat')
proj4string(pt1) <- CRS('+proj=longlat')
proj4string(pt2) <- CRS('+proj=longlat')
proj4string(pt3) <- CRS('+proj=longlat')
proj4string(pt4) <- CRS('+proj=longlat')
proj4string(pt5) <- CRS('+proj=longlat')
proj4string(pt6) <- CRS('+proj=longlat')
proj4string(pt7) <- CRS('+proj=longlat')
proj4string(pt8) <- CRS('+proj=longlat')
proj4string(pt9) <- CRS('+proj=longlat')
proj4string(pt10) <- CRS('+proj=longlat')
proj4string(pt11) <- CRS('+proj=longlat')

Web-mapping

library(maptools)
library(rgdal)
library(leaflet)
require(pals)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"), pt1$conductance_z)
m <- leaflet() %>%
  # Add tiles
  addProviderTiles("Esri.WorldTopoMap", group = "Topographical") %>%
  addProviderTiles("OpenStreetMap.Mapnik", group = "Road map") %>%
  addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
  addLegend(position = 'bottomright', opacity = 0.4, 
            colors = "blue", 
            labels = 'Participant 1',
            title = 'Cycling Routes (Netherlands)') %>%
  # Layers control
  addLayersControl(position = 'topright',
                   baseGroups = c("Topographical", "Road map", "Satellite"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  addCircles (data=pt1, group='Cycling routes', stroke = T, radius = 100, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt1$conductance_z),
              opacity = 0.2, color = conduct.pal(pt1$conductance_z)) %>%
  addLegend(values = pt1$conductance_z, pal = conduct.pal, 
            opacity = 1, title = "Skin Conductivity", position = "bottomleft")
m

Study Area

urban <-  c("Discontinuous dense urban fabric (S.L. : 50% -  80%)",  
            "Industrial, commercial, public, military and private units",
            "Discontinuous low density urban fabric (S.L. : 10% - 30%)", 
            "Isolated structures",
            "Continuous urban fabric (S.L. : > 80%)",
            "Discontinuous very low density urban fabric (S.L. : < 10%)",               
            "Discontinuous medium density urban fabric (S.L. : 30% - 50%)",  
            "Construction sites",  
            "Mineral extraction and dump sites")
natural <- c("Arable land (annual crops)",  
             "Pastures",   
             "Forests",    
             "Land without current use",
             "Green urban areas", 
             "Permanent crops (vineyards, fruit trees, olive groves)", 
             "Sports and leisure facilities",   
             "Herbaceous vegetation associations (natural grassland, moors...)",
             "Open spaces with little or no vegetation (beaches, dunes, bare rocks, glaciers)",
             "Wetlands")
roads  <- c("Other roads and associated land",
            "Railways and associated land",
            "Fast transit roads and associated land")
name = landcover$ITEM2012
landcover$group <- with(landcover, ifelse(name %in% urban, "urban",
                                          ifelse(name %in% natural, 
                                                 "natural", "roads")))
col = terrain.colors(3)
plot(landcover, col = col, col.regions = landcover$group,
                 edge.col = "transparent", axes = F,
                 colorkey = list(space = "bottom", height = 0.5, 
                                 width = 0.7),
                 main = "Study Area", sub = "Types of Landcover",
                par.settings = list(axis.line = list(col = 'transparent')))
legend("bottomleft", title = "Landcover", fill = col, 
       legend = c("Urban", "Natural", "Roads"))

Spatial Analysis & Visualization

Cycling Routes

# clipping landcover polygon to cycling route
neth_clipped <- as(crop(landcover, CHIPS_df), "SpatialPolygonsDataFrame")
col = terrain.colors(3)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"), pt1$conductance_z)
layout(matrix(1:6, nrow=1))
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt1, col = conduct.pal(pt1$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt2, col = conduct.pal(pt2$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt3, col = conduct.pal(pt3$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt4, col = conduct.pal(pt4$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt5, col = conduct.pal(pt5$conductance_z), cex = 1.5, pch = 20)
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt6, col = conduct.pal(pt6$conductance_z), cex = 1.5, pch = 20)

Space-time Cube

library(OpenStreetMap)
map <- openmap(as.numeric(c(max(pt1$lat), min(pt1$lon))),
               as.numeric(c(min(pt1$lat), max(pt1$lon))), type = "stamen-terrain")
transmap <- openproj(map, projection = proj4string(pt1))
map3d <- function(map, ...){
  if(length(map$tiles)!=1){stop("multiple tiles not implemented") }
  nx = map$tiles[[1]]$xres
  ny = map$tiles[[1]]$yres
  xmin = map$tiles[[1]]$bbox$p1[1]
  xmax = map$tiles[[1]]$bbox$p2[1]
  ymin = map$tiles[[1]]$bbox$p1[2]
  ymax = map$tiles[[1]]$bbox$p2[2]
  xc = seq(xmin,xmax,len=ny)
  yc = seq(ymin,ymax,len=nx)
  colours = matrix(map$tiles[[1]]$colorData,ny,nx)
  m = matrix(0,ny,nx)
  surface3d(xc,yc,m,col=colours, ...)
}
library(RColorBrewer)
# display.brewer.all()
bp = brewer.pal(11,"RdBu")
library(colourschemes)
cs = rampInterpolate(pt1$conductance_z, rev(bp))
pt1$conduct_pal <- cs(pt1$conductance_z)
plot3d(pt1$lon, pt1$lat, pt1$time, xlab="Longitude", 
       ylab="Latitude", zlab="Time", type = "s", 
       col = pt1$conduct_pal, size = 2.5, nam)

You must enable Javascript to view this page properly.